Spin nematic phases in models of correlated electron systems: a numerical study 
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Strongly interacting systems are known to often spontaneously develop exotic ground states under 
certain conditions. For instance, spin nematic phases have been discovered in various magnetic 
models. Such phases, which break spin symmetry but have no net local magnetization, have also 
been proposed by Nersesyan et al} in the context of electronic models. We introduce a A^- flavor 
microscopic model that interpolates from the large- A/" limit, where mean- field is valid and such a 
nematic phase occurs, to the more realistic = 1 case. By using a sign-free quantum Monte-Carlo, 
we show the existence of a spin nematic phase (analogous to a spin flux phase) for finite N] when N 
decreases, quantum fluctuations increase and this phase ultimately disappears in favor of an s-wave 
superconducting state. We also show that this nematic phase extends up to a finite critical charge 
doping. 

Dynamical studies allow us to clarify the Fermi surface property : in the nematic phase at half- 
filling, it consists of four points and the low-energy structure has a Dirac cone-like shape. Under 
doping, we observe clear signatures of Fermi pockets around these points. 

This is one of the few examples where numerical simulations show how quantum fluctuations can 
destroy a large- phase. 

PACS numbers: 71.27.+a71.10.Fd 



I. INTRODUCTION 

Spin nematic phases are characterized by an absence 
of local magnetic moment and long range orientational 
ordering of the spin degrees of freedom. Such phases 
were first introduced by Andreev and Grishchuk^i and 
then discovered in several magnetic modelsi^. Recently, 
spin nematic phases have also attracted attention in the 
the context of cold atoms where a gas bose condensate of 
^^Cr atoms in magnetic trap has been realized^. Due to 
the high spin, J = 3, of the ^^Cr atoms, it has been sug- 
gested that, upon release of the spin degrees of freedom, 
spin nematic states can be realized^'^. Mean- field real- 
ization of spin-nematic phases in particle-hole symmetric 
fermionic models on a square lattice have equally been 
proposed by Nersesyan et aim. Those phases are char- 
acterized by a checkerboard pattern of alternating spin 
currents around elementary plaquettes, and are coined 
spin flux phases (SFP). They break SU{2) spin symme- 
try as well as the lattice symmetry, but conserve time 
reversal symmetry. Such a phase could also possibly oc- 
cur in an extended Hubbard model^. The goal of this 
paper is to provide a microscopic realization of SFP. 

In this article, we concentrate on fermionic models on a 
square lattice and investigate the stability of the spin- flux 
phase. To do so, we consider the multi- flavored model 
Hamiltonian: 

H = 



(1) 




Here, the sum runs over nearest neighbors on a 
two-dimensional square lattice. The spinors c] ^ = 

(^1 T a' i a) ^^^^^ ^ ranges from 1 to N. In the limit 
A/" ^ oo, the saddle point approximation becomes ex- 
act, and we recover the mean- field results of Ref. 1. As 
N is reduced, quantum fluctuations around this saddle 
point are progressively taken in account and ultimately 
at A/" = 1 we recover a spin- 1/2 fermionic Hamiltonian. 
As we will see, the model, of Eq. (QJ allows sign free 
auxiliary field quantum Monte Carlo (QMC) simulations 
for arbitrary band fillings and values of N. Hence, we 
can map out the phase diagram as a function of doping 
for various A^, at fixed choice of coupling constant g/t. 
In the following, we fix t = 1 as the energy unit. 

The article is organized as follows. In the next section, 
we briefly describe the path integral formulation of the 
partition function which is the basis of the mean field 
and quantum Monte Carlo simulations. In section IIIII 
we present our numerical results which allow us to map 
out the phase diagram as a function of N and doping 
for a fixed coupling strength. Finally, we summarize our 
results in section HVI 



II. MEAN-FIELD AND QUANTUM MONTE 
CARLO SIMULATIONS. 

Both the saddle point equations and the QMC simu- 
lations are based on a path integral formulation of the 
partition function. Carrying out a Trotter breakup of 
the kinetic and interaction terms as well as a Hubbard- 
Stratonovitch transformation for the interaction term. 
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the partition function reads: 



I 



N 



where = ^cj^cj^^. Note that in the last equation, 

the trace runs over a single flavor. 

As is well known, in the large N limit the mean- field 
solution is recovered. With the substitution 



one obtains 



(3) 
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In this large N limit, the integral over the the fields rj which correspond precisely to the mean- field solutioni 
is dominated by the saddle point configuration and fluc- 



tuations around the saddle point are negligible. More- 
over, by neglecting the r dependence of the Hubbard- 
Stratonovitch fields, one obtains : 



A. Staggered spin flux mean fleld solution 

We consider the restricted order parameter: 



(7) 



(8) 



and the saddle-point equations read : 
dS{r^) 



= 



(5) This choice is certainly valid at half-band filling where 
perfect nesting pins the dominant instabilities to the 
wave vector Q = (7r,7r). Away from half- filling, order 
parameters at incommensurate wave vectors may be fa- 
vorable. This is not taken into account in the present 
Ansatz. With the above form, the mean- field Hamilto- 
nian reads: 



(6) 
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where the sum over k is restricted to the magnetic Bril- 
louin zone (MBZ) and with 

z{k) = -2i (sin(A:^ + Q^/2) - sin(A:^ + Qy/2)) (10) 
Thus, we obtain the saddle point equation 

^= ^ Yl ^(4,a'4+Q,j 
keMBZ,a 

( ^f](/J': )) (11) 

y az{k) J \ ^k+Q,c7 J 
which we solve self-consistently. 

Mean Field 
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FIG. 1: (Color online) Mean Field order parameter as a func- 
tion of band-filling for various coupling strengths g. 

Figure Q shows the order parameter as a function of 
doping. The staggered spin flux phase survives up to a 
finite critical doping, where a first order transition to the 
paramagnetic phase occurs. This transition is signaled 
by a jump in the order parameter. At half-band filling, 
the single particle dispersion relation is given by, 

E{k) = ±^e^{k)^A^{k) (12) 

with A(/c) = 2gr]{cos{kx) — cos{ky)). Hence, it exhibits 
Dirac cones around the (±7r/2, ±7r/2) /c-points and the 
Fermi surface is given by those four points. In the dop- 
ing range where the order parameter does not vanish the 
Fermi surface consists of hole-pockets centered around 
the above mentioned ^-points (See Fig. Note that 
due to the d-wave symmetry of the order parameter, the 
Fermi surface is not invariant under reflections across the 
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(1,1) axis (i.e. kx ky and ky kx). As appropri- 
ate for a first order transition, the Fermi surface changes 
abruptly from hole pockets around (7r/2, 7r/2) in the spin 
flux phase to a large Fermi surface centered around (0, 0) 
in the paramagnetic phase. 



B. Quantum Monte-Carlo 

The Hamiltonian occurring in is quadratic in the 
fermionic variables so that the trace can be carried out 
analytically to obtain: 

j D^e-^i^^o),r^l.j)^^y^[deiM{^)f (13) 

Since the spin current is even under time reversal sym- 
metry, it has been shown that, when ^ > 0, the fermionic 
determinant is positive for each Hubbard-Stratonovitch 
configuration^. This absence of sign problem is valid for 
any lattice topology and any filling. Hence each config- 
uration of Hubbard Stratonovitch fields can be sampled 
according to its weight with Monte Carlo techniques. 

In the following, we study the ground-state proper- 
ties of the Hamiltonian with the projector auxiliary 
field QMC algorithm on a two-dimensional square lattice. 
The details on this algorithm can be found for example in 
Ref. 10, where the authors consider a very similar model 
from the technical point of view. Dynamical informa- 
tion is obtained by using a recent implementation of the 
stochastic analytical continuationiii^. 



III. NUMERICAL RESULTS 

A. Phase diagram 

For = 1, the interaction term of our model can be 
rewritten, up to a constant as: 




= - 3^^(4^4^c^.^c^.^ + /i.c.) 

- 2gY,S,-S, (14) 



4 



n(k), g=.1,n = 0.74 



functions: 
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FIG. 2: (Color online) Mean-field solution of n{k) over the 
Brillouin zone as a function of band filling for a fixed g — 1. 



As apparent, there is a pair- hopping term so that we 
can expect an on-site superconducting instability which 
we pick up by measuring the equal time pair correlation 



sciR) = {cyyoiCo,}- 



(15) 



In the presence of long range off-diagonal order, its 
Fourier transform at Q = (0, 0) 



5C(Q = (0,0)) = -1^5C(^) 



(16) 



should converge to a finite value in the thermodynamic 
limit, whereas it vanishes as in the absence of long- 
range order. 
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FIG. 3: (Color online) Finite-size scaling of the supercon- 
ducting correlations vs the inverse number of sites of 
the system (L is the linear size) for the half-filled case and 
different N. The Q = (0,0) Fourier correlations are plotted 
(open symbols) and, when long range order is present, we 
also show the largest distance real-space correlations (filled 
symbols) . 

The order parameter of the spin flux phase expected 
in the limit of A/" ^ oo, reads: 



J^(0,0 + x)- J^(0,0 + ^)1) 



where 



r{R,R^x) 



it,/i ^ R-\-x^v 



(17) 



(18) 



is the bond spin current operator along x (summation 
over repeated indices is assumed). A similar defini- 
tion holds along the y direction. As expected from the 
mean-field solution, we have observed that the strongest 
signal occurs in the d-wave symmetry channel and the 
Fourier transform of these correlations is maximum at 
Q (7r,7r). 

Technical details: we are able to simulate up to 20 x 
20 square lattices with a projection parameter 9t = 10. 
We fix Ar = 0.1 for the Trotter decomposition. For 



5 




0-ON=1 
□-□N=2 
<^ N=3 

A^N=3 (L/2,L/2) 
<^ N=4 

N=4 (L/2,L/2) 



0.01 0.02 0.03 0.04 0.05 0.06 0.07 

FIG. 4: (Color online) Similar as the previous plot for d-wave 
spin current correlations. 
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FIG. 5: Order parameter 77 as a function of N at half- filling. 
The data point dX N — 00 corresponds to the mean-field value. 
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FIG. 6: (Color online) Finite-size scaling of the spin-spin cor- 
relations vs the inverse volume, l/L^, of the system for the 
half-filled case and different N. For all cases, the Q — (7r,7r) 
Fourier correlations (which is the largest one) vanish in the 
thermodynamic limit, indicating the absence of magnetic or- 
der. 



simplicity, in the following, we fix the coupling constant 
5=1- 

Half-filling case : On Fig. [31 we plot the scaling of the 
superconducting correlations for different values of TV. 
According to our definition ([TB|) . the Q = (0,0) correla- 
tions are divided by the number of sites so that a finite 
value in the thermodynamic limit indicates long range 
order. Moreover, when this is the case, we also plot the 
largest distance real-space correlations that should con- 
verge to the same value in the thermodynamic limit. In 
particular, we have a clear signal of s-wave superconduct- 
ing phase for A/" = 1 and N = 2 but this phase disappears 
for larger N . 

At = 1 and under the canonical transformation, 

4^(-irc,,i 4^ ^(-1^4^ (19) 

the interaction term transforms as, 



-|(n, - l)(n, - 1) +^(4T4,iS-,iS-,T + ^•^•)' (^0) 

and the kinetic energy remains invariant. The observed 
long-ranged pairing correlations map onto long ranged 
tranverse spin-spin correlations. Since the transformed 
model has an SU (2) spin symmetry, the transverse spin- 
spin correlations on any finite lattice take the same value 
as the longitudinal spin-spin correlations. Transforming 
back to the original model maps the longitudinal spin- 
spin correlations to a charge density wave order. Hence 
at half-filling and A" = 1 pairing correlations and charge 
density wave correlations are locked in by symmetry. 
Clearly, doping breaks this symmetry^^. 

Similarly, we investigate the occurrence of SFP 
phase by computing the corresponding correlations (see 
Eq. (|T7|) ). Moreover, for a finite system, the SU{2) sym- 
metry of these correlations cannot be broken so that we 
average over the 3 components to reduce the error bars. 
For example, we plot on Fig. (b) all three components 
of the SFP correlations showing that indeed SU (2) sym- 
metry is restored (within the error bars) despite being ex- 
plicitely broken for a given Hubbard- Stratonovitch con- 
figuration. On Fig.^ we plot the scaling of these correla- 
tions showing that SFP is present in the thermodynamic 
limit for A^ > 3 but is not stable for smaller values of A'. 
Again, when the Q = (tt, tt) Fourier component indicates 
long range order (LRO), we also plot the scaling of the 
largest distance correlations that should have the same 
value in the thermodynamic limit for an ordered phase. 
In particular, we recover the existence of SFP phase for 
large A" as found at the mean-field level. A non trivial 
result is that this phase survives for finite N > 3 since 
our QMC simulations, which are free of the sign prob- 
lem, include quantum fluctuations. We observe rather 
strong finite size effects and in particular, the signal is 
weaker when the lattice contains (7r/2,7r/2) k-points in 
its Brillouin zone. Indeed, as we will discuss in the follow- 
ing, this point correspond to the low-energy excitations. 
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FIG. 7: (Color online) (a) Spin current correlations, averaged 
over the 3 components, vs k along the (TXM) path in the 
Brillouin zone. Parameters are L = 12, iV = 10 and g = 1. 
As a function of doping 6 (increasing from top to bottom), 
the maximum at (tt, tt) is rapidly suppressed, (b) Same for 
L — 10 and N — 3. Here, we plot all 3 components showing 
that SU{2) symmetry is recovered in our data. As a function 
of doping (increasing from top to bottom), the maximum at 
(tt, tt) is rapidly suppressed. 



Since the staggered spin current correlation function of 
Fig. ^converges to the square of the order parameter 77, 
we can extract this quantity as a function of N. Our re- 
sults are plotted in Fig.O As apparent, the Monte Carlo 
results at finite values of N smoothly scale to the mean 
field result valid at A/" = 00. 

In the SFP at A" > 3, we have computed the spin- 
spin correlation functions (5^- Sj). As shown on Fig. El 
our results show no sign of long-range spin ordering and 
hence the absence of a magnetic moment. This confirms 
the point of view that the spin- flux phase that we observe 
for A^ > 3 indeed corresponds to a spin nematic phase. 

Doping : Since the SFP instability occurs in the 
particle-hole channel, it is favored at half-filling where 
the Fermi surface exhibits perfect nesting. But the ques- 
tion about its existence under doping remains open. Our 
QMC simulation being free of the sign-problem for any 
filling, we have computed the phase diagram as a func- 
tion of doping for different values of A'. As expected 
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FIG. 8: (Color online) Phase diagram as a function of doping 
and 1/N obtained with QMC simulations. At large N, we 
recover the existence of a spin flux phase (SFP) over a finite 
doping region. For larger doping or smaller N, the system 
develops s-wave superconducting correlations (SC) or stays 
in a metallic phase (see text). 



from the mean-field results, we show on Fig. [HI that the 
SFP phase has a finite extension up to a critical doping. 
In particular, when N = 10, we find a critical doping 
around 0.15. 

On Fig. 13 we plot the spin current correlations vs k 
along the (TXM) path in the Brillouin zone (shown in 
Fig.[H|). For doping S smaller than a critical value, we see 
a large maximum at (tt, tt), whereas for larger doping, the 
correlations become featureless, indicating the disappear- 
ance of SFP order. We also observe that the maximum of 
the spin-current correlations move away from (tt, tt): for 
example, when L = 12, A^ = 10 and with doping ~ 0.17, 
the maximum is at (57r/6,7r) and equivalent k-points. 

By performing a finite-size scaling analysis of our data, 
we are able to draw conclusion about the existence or not 
of a SFP in the thermodynamic limit as a function of dop- 
ing and A". These results are summarized in the phase 
diagram shown in Fig. [HI Several comments are in order : 

i) Under doping the SFP only survives for large values of 
A", hence showing that it is very sensitive to fluctuations. 
At A" = 10 the critical doping, 5c — 0.15 at which the 
SFP vanishes is substantially smaller than the N ^ 00 
restricted mean- field result Sc — 0.25. In the vicinity of 
the phase transition from the SFP to the paramagnetic 
phase transition at A^ ^ 00, fluctuations will play an 
important role. Hence at large but finite values of A', 
substantial differences can be expected. Furthermore, 
our restricted mean-field solution does not allow for in- 
commensurate ordering as suggested by the Monte Carlo 
results. 

ii) The rest of the phase diagram is dominated by an 
s-wave SC phase. At A" = 1, the s-wave SC persists of 
course away from half-filling since it is a particle-particle 
instability that does not require any nesting property. As 
A" grows, the s-wave SC order parameter is reduced and 
ultimately vanishes to produce the paramagnetic phase 
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in the AT ^ oo limit. At = 10 and 5 > 5c ^ 0.15 
we observe no SFP correlations and no SC correlations 
either. Therefore, it seems that we might have a metallic 
phase with a Fermi surface. However, due to the presence 
of pair hopping in the model, we know that this metallic 
phase will ultimately be unstable at very low tempera- 
ture toward s-wave superconductivity. This instability 
probably occurs at too low temperature to be observed 
in our simulations. To confirm this point of view, let us 
write the interaction term as: 



N 




b (3=1 



(21) 



where b = (z, j) denotes a nearest neighbor bond, and 
J^Q, is the spin current on bond b for flavor index a. This 
forms shows that in the large-N limit, the flavor index /3 
is embeded in the mean field spin current of the other fla- 
vors, ^ T^c^^p ^5%. The second term, j^J^^p-J^^p is a fac- 
tor 1/N smaller than the mean- field term and provides 
the pair hopping term as explicitly shown in Eq. ([T^ . 
At = oo, only the mean- field term survives. For large 
but finite values of A^, a small pairing term is present and 
will in the paramagnetic phase trigger a superconducting 
instability at low temperatures. 

iii) The coexistence a spin flux phase and superconduc- 
tivity at finite doping cannot be excluded. As we will 
see below, doping the SFP leads to a Fermi surface con- 
sisting of hole pockets centered around the (±7r/2, ±7r/2) 
points in the Brillouin zone. Due to the presence of pair 
hopping, this Fermi surface should be unstable towards 
a superconducting state. As a result, for intermediate 
values of A^ such that we observe SFP at zero doping 
and such that the superconducting signal at finite dop- 
ing is not strongly reduced due to a 1/N factor, we have 
found numerical evidence that supports coexistence of 
both phases for small doping (for instance, A^ = 4 and 
S 0.06 in the phase diagram). 



B. Dynamical properties 

To study the charge degrees of freedom, we compute 
the single particle spectral function A{k^ uj) which is re- 
lated to the imaginary time Green function via: 



1 

(ci (r)cr: ) = - / 



dcje-^^A(/c, - 



(22) 



In order to plot this single particle spectral function, 
we follow a path in the Brillouin zone shown on Fig. El In 
Fig. El we first plot two examples of A{k^u;) in the SFP 
at half-band filling. As expected from the mean-field cal- 
culation we observe a clear sign of the Dirac cones around 




FIG. 9: In the following spectral function, from bottom to 
top, the k value follows the path F, X, M and then, along 
the (tt, 0) to (0, tt) line. 
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FIG. 10: (Color online) Spectral functions vs frequency ou 
for different k points at half-filling. Up: N = 10 and L — 
12; bottom: N — A and L — 16. Both set of parameters 
correspond to SFP and Dirac cones around (7r/2,7r/2) are 
clearly seen. 



the (7r/2,7r/2) /c-points in the Brillouin zone. The overall 
dispersion relation compares favorably to the mean-field 
result 



E{k) = ±Je^{k)^A^{k) 



(23) 



with A(^) = 2gr] {cos{kx) — cos{ky)). 

In order to be more quantitative, we plot the dispersion 
relation of the main branch (by looking at the maximum 
of the spectral function) as a function of the distance to 
the nodal point. On Fig. ^2 we clearly see the Dirac cone 
structure with different velocities parallel and perpendic- 
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FIG. 11: (Color online) Dispersion relation obtained with 
N = 4 and L = 16 at half-filling. The energy is taken at 
the maximum of the spectral function and data are plotted 
vs distance to (7r/2, 7r/2). The small asymmetry is due to the 
uncertainty in the spectral functions obtained with a Maxi- 
mum Entropy technique. 
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FIG. 12: (Color online) Spectral functions vs a; — /x for dif- 
ferent k points. Parameters are: = 10, L = 12 and doping 
0.12. It corresponds to a doped SFP phase and Dirac cones 
around (7r/2,7r/2) are clearly seen. The chemical potential is 
= -0.48. 



ular to the Fermi surface. By fitting the slopes, we ob- 
tain for this set of parameters, v± =3.11 and v\\ = 0.578. 
At the mean- field level, according to Eq. (|23j) . these two 
velocities are equal respectively to 2^/2 and 2\^gr] so 
that their ratio gives directly access to the SFP order 
parameter r] (when g = 1). Of course, the QMC data 
include some renormalization of these values and the ra- 
tio of v\\/v± = 0.186 is inconsistent with the SFP order 
parameter ( rj 0.45) as extracted from the spin current 
correlations (see Fig.0}. 

In the mean-field approach, doping the SFP leads to 
hole pockets centered around the (7r/2,7r/2) ^-points. 
Numerical QMC simulations confirms this as shown on 
Fig. El At A/" = 10 and S = 0.12 we observe a clear 
signature of the Dirac cone, however the quasiparticle at 
k = (7r/2,7r/2) lies above the Fermi energy. 
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FIG. 13: (Color online) Spectral functions vs frequency cu for 
different k points. Top: N = 3 and L = 12 for doping 0.12 
(/i = -0.7). Bottom AT = 2 and L = 16 at half- filling (/i = 0). 
Both set of parameters correspond to s-wave superconducting 
phase. 



Deep in the superconducting phase, the single particle 
dispersion relation follows the the mean-field form: 



E{k) = ±V(e(fc)-Ai)2+A2, 



(24) 



with /c-independent superconducting gap. Such a disper- 
sion relation is consistent with the Monte Carlo data of 
Fig. [iniat = 3 and S = 0.12. In the proximity of 
the SFP at N = 2 and half-band filling (see Fig. ^ the 
dispersion is not captured by the above form. In particu- 
lar, a modulation of the gap along the (tt, 0) to (0, tt) line 
is observed. This modulation is reminiscent of the Dirac 
cone structure observed in the spin- flux phase, and hence 
allows the interpretation that short range spin currents 
survive in the superconducting state in the proximity of 
the SFP. 



IV. CONCLUSION 

We have shown the existence of a spin nematic phase 
in a two-dimensional electronic model. This phase had 
been proposed at the mean-field level and it is the aim 
of this article to investigate its stability against quan- 
tum fluctuations. In order to do so, we have performed 
unbiased QMC simulations (without sign problem) for a 
variety of models with a flavor parameter N. With this 
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trick, we can interpolate from the large N limit where the 
mean- field is valid to finite N regime where quantum fluc- 
tuations around the saddle point are progressively taken 
into account. We find that the spin- flux phase is stable 
for a range of finite N but is ultimately destroyed for the 
more realistic = 1 case. This is one of the few examples 
where QMC simulations are able to show how quantum 
fluctuations can destroy a large- A/" phase^^. Since our 
model includes pair hopping processes, when the spin- 
flux phase is found to be unstable, it is replaced by an 
on-site 5-wave superconductivity. 

We have also investigated the phase diagram as a func- 
tion of doping since our QMC simulations are free of the 
sign problem for any filling. Because the spin-flux phase 
is due to the nesting property of the Fermi surface, it is 
expected to disappear with doping. This is indeed what 
we found but still, this phase can persist over a finite 
range of doping, similarly to what is found at the mean- 
field level. 



Using a recently developed analytical continuation 
technique, we have computed dynamical spectral func- 
tions. In the SFP phase, we clearly see Dirac cones form- 
ing around (7r/2, 7r/2) and equivalent /c-points, where the 
dispersion becomes hnear. Up to a critical doping, these 
structures are stable so that the Fermi surface becomes 
pocket- like. 
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